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We consider perhaps the simplest, nontrivial problem in advection-diffusion — a finite 
absorber of arbitrary cross section in a steady two-dimensional potential flow of concen- 
trated fluid. This problem has been studied extensively in the theory of solidification 
from a flowing melt, and it also arises in advection-diffusion-limited aggregation. In both 
cases, the fundamental object is the flux to a circular disk, obtained by conformal map- 
ping from more complicated shapes. Here, we construct an accurate numerical solution 
by an efficient method that involves mapping to the interior of the disk and using a 
spectral method in polar coordinates. The method combines exact asymptotics and an 
adaptive mesh to handle boundary layers. Starting from a well-known integral equation 
in streamline coordinates, we also derive high-order asymptotic expansions for high and 
low Peclet numbers (Pe). Remarkably, the "high" Pe expansion remains accurate even 
for such low Pe as 10^"^. The two expansions overlap well near Pe = 0.1, allowing the 
construction of an analytical connection formula that is uniformly accurate for all Pe and 
angles on the disk with a maximum relative error of 1.75%. We also obtain an analytical 
formula for the Nusselt number (Nu) as a function of Pe with a maximum relative error of 
0.53% for all possible geometries after conformal mapping. Considering the concentration 
disturbance around a disk, we find that the crossover from a diffusive cloud (at low Pe) 
to an advective wake (at high Pe) occurs at Pe « 60. 



1. Introduction 

The transfer of mass, heat, or other passive scalars in fluid flows is a major theme in 
transport science (Lea l, 1992). The canonical model problem involves a uniform back- 
ground flow of speed, C/qo, and concentration. Coo, past an absorbing object of charac- 
teristic size, L. Given the steady, incompressible flow field, u (scaled to Uoo), the steady 
tracer concentration around the object, c (scaled to Coo), satisfies the (dimensionless) 
linear advection-diffusion equation, 

PeoM-Vc = V^c, (1.1) 

where Peo = UooL/ D is the bare Peclet number, which measures the relative importance 
of advection compared to diffusion (with a diffusivity D) . The partial differential equation 

f The computer programming codes used for the results in this paper are available at 
http : //www . advection-diffusion. net. 
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C = U — X, |z| — > (X) 



Figure 1. The fundamental problem of advection-diffusion; A finite absorber (source) in a 
steady flow of uniformly concentrated (depleted) fluid in the complex z plane; fiz is the region 
exterior to the object. 

(PDE) Hl.l() must be solved subject to the boundary conditions (BCs) c = on the object 
and c = 1 far away (at oo), to obtain the dimensionless normal flux density, 

a^fi-Vc, (1.2) 

everywhere on the surface of the object. (Since the problem is linear, we may also consider 
the equivalent problem of a source object at c = 1 relative to a depleted background fluid 
at c = 0.) 

Although similarity solutions exist for infinite leading edges (see below), in the usual 
case of a finite absorber, the mathematical problem is intractable. It can even be difficult 
to solve numerically because Peo appears as a singular perturbation in the PDE in both 
limits, Pco — > oo and Pco 0. The classical approach, therefore, has been to employ 
asymptotic analysis to obtain approximate solutions, usually relating the total integrated 
flux, or Nusselt number, Nu, to the Peclet number, Pe. Early studies of this type focused 
on spheres ( |Acrivos fc Taylor! 1962) and more complicated shapes (IBrennerl 1963) in 
Stokes flows, and later studies dealt with heat or mass transfer in steady shear flow 
| |Phillips[ 1990). Related work continues to the present day, e.g. in the context of nutrient 
uptake by single-cell organisms ( |Magar et al.', 2003). 

The Nu(Pe) relation contains useful global information, but one sometimes requires a 
complete solution to the problem, including the local flux profile on the absorber. In this 
paper, we focus on a well known special case, ideally suited for mathematical analysis and 
physical interpretation: steady advection-diffusion in a two-dimensional, irrotational flow. 
The velocity field is described by the flow potential, 0, as it = V(/) and the dimensionless 
boundary-value problem (BVP) is: 

Pe„V(/. • Vc = V\ V'^-O, {x,y) € VL, (1.3) 
c=l, n-V(/) = 0, {x,y) e dVlz and c -> 0, -f oo; (1.4) 
V(/)->&, + ^ oo, (1.5) 

where 51^ is the flow region, exterior to the object's boundary, d^z, as shown in Fig. ^ 
Note that in 11. 4|) we use boundary conditions for desorption (c = 1 on the object and 
c = far away), which are somewhat more convenient that those of adsorption (c = 
on the object and c = 1 far away). The two problems are mathematically equivalent by 
linearity. 

The key to analyzing (|1.3|l is to view points in the plane as complex numbers, z = 
x + iy. For example, this enables a transformation to streamline coordinates, $ = (ji + iip, 
which reduces the problem to a thin absorbing strip in a uniform flow ( |Boussinesq| 
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1905). For a more general perspective in terms of conformal mapping, see Bazant (2004). 
Using complex-variable techniques, the general system of PDEs H1.3|l has been studied 
recently with applications to tracer dispersion | |Koplik et al.\ 1994) and heat transfer 

(Morcga & Bejan' 1994) in porous media, as well as vorticity diffusion in strained wakes 
Eames & Bush, 1999: IHunt fc EamesI 2002). 

Complex analysis becomes particularly useful when the interface, dflz(t), is a mov- 
ing free boundary, driven by the local flux density, CTz. For a broad class of transport- 
limited growth phenomena, the interfacial dynamics, whether deterministic or stochas- 
tic, can be formulated in terms of a time-dependent conformal map from a simple, 
static domain to the evolving, physical domain IjBazant et aL\ 2003). For continuous 
growth by advection-diffusion, this approach was introduced by |Wijngaarden| (1966) and 
IMaksimovl ('1976'). who solved H1.3(l - (|1.5(l in streamline coordinates. Maksimov's method 
has been used extensively by Kornev and his collaborators to describe solidification and 

freezing from a flowing melt fChugunov fc KornevI 1986; 'Kornev k Chugun ov| 1988; 

IKornev fc Mu khamad ullinal 1994: .Alimov ef"oi.l 1994: |Alimov et al.i 1998; |Cummings fc Kornev| 
1999; [Uummin gs et aZ.J 1999). These studies significantly extend the conformal-map dy- 
namics for Laplacian growth by pure diffusion, without advection IjPolubarinova-Kochinal 
1945: l(;alinl 1945: iHo^iscml 1992). 

Here, we are motivated by a new, discrete growth model, Advection-Diffusion-Liniited 
Aggregation ( ADLA) , which describes the growth of fractal aggregates in a fluid flow via 
a stochastic conformal map (|Bazant et al.\ 2003). In that case, the BVP (|1.3(l - p.5|l must 
be solved for a circular absorber for all values of Pcq. The normal flux distribution, a^, 
determines the probability measure for growth events. The detailed description of the 
growth measure in this paper shows how dynamics crosses over from diffusion-dominated 
to advection-dominated growth regimes, as a function of the time-dependent Peclet num- 
ber. These two "fixed points" of the dynamics are related to special similarity solutions 
l|Bazantl 2004), which correspond to the asymptotic regimes of high and low Pco analyzed 
in this paper. 

Beyond such applications, the BVP l|1.3|l - l|1.5|l . for an arbitrary single-connected do- 
main, merits serious mathematical study in its own right. It is perhaps the simplest 
advection-diffusion problem with a nontrivial dependence on the Peclet number. It may 
also be the most complicated advection-diffusion problem for which a nearly exact ana- 
lytical solution is possible, as we show here. 

The paper is organized as follows. In Sec. |2] we set the stage for our analysis by 
reviewing two key properties of the BVP H1.3(l - H1.5|l : (i) It can be recast as a singular 
integral equation in streamline coordinates, and (ii) conformal mapping can be applied to 
work in other convenient geometries for numerical solution and mathematical analysis. In 
Sec. 121 we present an efficient, new numerical method to solve the BVP, after conformal 
mapping to the interior of a circular disk. In Sec. 01 we derive accurate asymptotic 
expansions for a for high Peclet numbers by applying an exact iterative procedure to the 
integral equation. In SecElapproximate formulae for a are derived when Pe is sufficiently 
small via approximating the kernel of the integral equation and solving the resulting 
equation exactly by known methods. In Sec. Elan accurate ad hoc connection formula is 
given for a by combining the formulae for high and low Pe, and it is also integrated to 
obtain the Nu(Pe) relation. Finally, in Sec. [3 we conclude with a discussion of some of 
the implications and applications of our results, as well as by posing a few challenges for 
future work. 
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2. Mathematical Preliminaries 

2.1. Streamline Coordinates 

In his analysis of high Reynolds number flows, |Boussinesq| (1905) discovered a hodo- 
graph transformation (exchanging dependent and independent variables) which greatly 
simplifies (|1.3|l . It is well known that the velocity potential is the real part of an analytic 
complex potential, $ = + if/;, where the harmonic conjugate, is the stream func- 
tion IjBatchelorl 1967) . Using the Cauchy-Riemann equations, it is easy to show that the 
concentration profile, c, satisfies the simplified PDE 

_ ^ 

° d(t> dcj)^ dip"^ ' 



^ dc d c d c 

PeoTTT = TTT? + (2-1) 



in "streamline coordinates", (0, V')- The physical significance of this equation is that 
advection (the left hand side) only occurs along streamlines, while diffusion (the right 
hand side) also occurs in the transverse direction, along isopotential lines. 

Boussinesq's transformation corresponds to a conformal mapping to a plane of a uni- 
form flow, described by a constant V0. Therefore, an arbitrary domain, fi^, as shown 
in Fig. n is mapped to the exterior of a straight line segment, or strip, parallel to the 
streamlines (which is a branch cut of the inverse map). Some examples are given in Fig. El 
In streamline coordinates, the BCs (|1.4|) and H1.5|l are transformed as follows: 

c=l, tp ^ 0, -2A < (j)< 2A, and c ^ as 0^ + -^^ ^ oo. (2.2) 

The constant A is determined so that 4A is equal to the difference of the flow potential 
between two stagnant points, for example, the points Su and Sd shown in Fig. ^ 

2.2. Formulation as an Integral Equation 

In streamline coordinates, the advection-diffusion process past a finite strip can be for- 
mulated in terms of an integral equation using the classical method of Green's functions 
( |Stakgold, 1998). For reasons to become clear in Sec. 12.31 we let a; = 0, j/ = V', ^ = 1/2, 
and Pe = Peo/2, so the EVP (gUJ-lEiSl takes the form 

dx dx^ ^ dy-^ ' ^ ^ 

c = 1, on y = 0, — 1 < a; < 1 and c = 0, as -I- — > oo. (2-4) 

Green's function G{x,y) for (|2.3|l and 1)2.4(1 . which expresses the concentration profile 
generated by a unit source flux at the origin, satisfies the PDE 

2Ve—- — - — = 5[x)5{y) 
dx dx"^ dy^ 

After removing the first derivative term by a change of variables and using polar coordi- 
nates, 

G(a;,y) = eP^-F(r,0), (2.5) 
we find that F obeys the Helmholtz equation 

- Pcp]F = 5{x)5{y) , (2.6) 

whose solution is a modified Bessel function of the second kind, K^^Per), where r = 
(a;^ +y^)^^^- Taking into account the unit normalization, we obtain Green's function, 

G(x,2/) = eP^"i^o(Per)/2^. (2.7) 
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Figure 2. Numerical solutions using the method of Sec. El for the concentration profile (color 
contour plot) and streamlines (open yellow curves) around different absorbing objects in a 
uniform background potential flow: (a) the unit disk, (b) a finite strip, (c) a square, and (d) an 
ADLA fractal cluster with Peo = 1, 2, 1.2 and 0.05 respectively. The geometries in (b), (c) and 
(d) are obtained by conformal mapping from (a). Case (b) corresponds to streamline coordinates; 
case (c) is obtained by the numerical Schwartz-Christoffel mapping (,Trefethent 1986); case (d) 
is obtained from a stochastic, iterated conformal map (Bazant et aQ 2003). The horizontal axis 
is labeled by numerical values of the spatial coordinate x. However, the distances in the figures 
(a)-(d) are scaled by the "conformal radius" Ai, see H2.11|l below, which is the characteristic size 
of the physical object. The Ai is chosen so that the renormalized Peclet number is the same, 
Pe = ^iPeo — 1, in all cases, which explains why the far-field solutions look the same. 



The concentration profile, c, everywhere in streamline coordinates is obtained by con- 
volving Green's function, G, with the flux on the strip, 

c{x,y)^J G{x-x',y)-2a{x')dx', (2.8) 

where the factor 2 is included because the flux has the same value, 

,(,)^_|£(,,o+)^|(x,o-), 

on the upper and the lower sides of the strip, respectively. Therefore, the boundary value 
problem described by (|2.3|l and l|2.4|l is equivalent to finding the (t{x) that satisfies the 
integral equation ^Pearson. 1957; , Wij ngaarden[ 1966) 

J e^''^''-'^">Ko{Pe\x-x'\)a{x')dx' = TT, -1< x < 1, (2.9) 
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which forms the basis for the theory of soUdification in flowing melts (|Maksimovl, 1976). 
In this context, H2.9|l has been analyzed for large and small Pe by Kornev and collabora- 
tors, as cited in the introduction. Below, we will extend these results and construct an 
analytical approximation that is uniformly accurate in both Pe and x. 

The reader may worry about the existence and uniqueness of solutions because l|2.9|l is 
a Frcdholm-type equation of the first kind (with a difference kernel). In the present case, 
however, the symmetrized kernel, Kq{Fc\x ~ x'\), is positive definite, and thus invertible; 
see Appendix IXI for an explanation. 

It seems tempting to approach (|2.9|l using Fourier-type methods because it involves 
a convolution IjTitchm arsh 1948), but the kernel is not a periodic function. Wc also 
note that the kernel is singular and not of the classical Cauchy type ((Muskhe lishvilii 
1992). It is known that l|2.9|l admits a solution which can be expanded in terms of Math- 
ieu functions IjRvachev ' 1956:' Protsenko fc Rvachevl 1976), but such representations are 
impractical for computations over a wide range of Pe, and give no insight into the de- 
pendence of the fiux cr on a; and Pe. The collocation method has been used successfully 
to obtain numerically the solution of H2.9|l IjKornev fc Mukhamadullinal 1994). 



2.3. The General Principle of Conformal Invariance 

There is a simple way to understand why Boussinesq's transformation works: The advection- 
diffusion PDE (|1.3() is invariant under conformal changes of variables, even though its 
solutions are not harmonic functions, which also holds more generally for some other 
equations IjBazantl 2004). As such, the boundary value can be transformed to any con- 
venient geometry by conformal mapping. Streamline coordinates is a good choice for 
asymptotics, but other choices are better suited for numerical analysis and similarity 
solutions. 

Here, we exploit this general principle to map the BVP H1.3|l - H1.5|l to other useful 
coordinate systems. If — Kc^{w) and — F{'w,W) (where W denotes the complex 
conjugate of w) so lve (jl.^ll in some simple domain, fi^j, then (f>z — Re$(/(z)) and 
Cz = F{f{z),f{z)) solve ()1.3|l in an arbitrary mapped domain, U.^ = g{U,z), where 
z = g{w) — f^^{w). The concentration BCs (|1.4|l are conformally invariant, but, since 
the BC (|1.5ll prescribing the background fiow is not, due to the gradient, care must be 
taken in transforming the solution. 

For advection-diffusion-limited growth (IBazant et al.\ 2003) , it is natural to let fi^ be 
the exterior of the unit circle, so our canonical problem is that of a concentrated flow 
past a circular absorber, shown in Fig.|3a), with the velocity potential 

(/) = Re + |w| > 1. (2.10) 

To reach other geometries, the mapping, z — g{w), must be univalent (conformal and 
one-to-one), so it has a Laurent series of the form 

g{w)^Aiw + Ao + — + ..., \w\>l, (2.11) 
w 

where Ai is a positive, real constant which defines an effective diameter of (the 
"conformal radius"). In order to preserve the BC H1.5|l which sets a unit fiow speed at 
z = oo, we would need to set the dimensionless fiow speed to Ai at w — oo. Instead, 
we choose to redefine the velocity potential in riu, to preserve the unit flow speed, as in 
H2.10|l . and then define a renormalized Peclet number, Pe = Ai Peg, analogous to the 
time-dependent Peclet number for ADLA defined bv lBazant et aLI ('2003'). 
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The BVP for the concentration flw then becomes 

PeV0-Vc=V2c, |w| > 1, (2.12) 

c = 1 on = 1, and c — as luij ^ oo. (2-13) 



The physical significance of the renormalized Peclet number, Pe, is that it determines 
the far-field solution, independent of the absorber's shape and the bare Peclet number. 
This point is illustrated in Fig. |21 where the concentration and flow field far away from 
various objects at Pe = 1 looks the same, in spite of extremely different shapes, ranging 
from a circle to a fractal ADLA cluster. Therefore, we view Pe as the basic parameter in 
our analysis, from which we define the bare Peclet number, Pcq = Pe/Ai, for arbitrary 
domains, in terms of the univalent map from the exterior of the unit circle. 

From this perspective, streamline coordinates are obtained via the Joukowski transfor- 
mation IjCarrier et al.\ 1983), g{w) — {w + l/w)/2, which maps the unit circular disk onto 
the finite strip of length 2 centered at the origin along the real axis. The BVP (|2.12l) - 
(|2.13|l is then transformed to the form H2.3|l - I|2.4|l given above. In general, the fluxes on 
the boundaries dilz (the absorber surface) and dftw (the unit circle, w ~ e'^) are related 

by 

a,„(0;Pe) = \gfiw)\ a,ig{w);Fe), (2.14) 

where is the flux on dflw and cr^ is the flux on dflz- In the case of streamline coordi- 
nates (|2. 311 - 1)2. 4|l . the flux on the strip, az{x; Pe), is thus related to the flux on the circle, 
a»(0;Pe), by 

cr^(6»;Pe) = |sin6l| ct^(cos6I; Pe). (2.15) 

For a bounded flux on the circle, aw, the flux on the strip, az, always diverges as 0[{1 — 
3;2^-i/2j ^ approaches ±1. Therefore, although we will use the strip geometry, Qz, for 
asymptotic analysis in Sees. 0] and [SJ the circle geometry, 57^,, is a better starting point 
for our numerical analysis in Sec.|3| In all cases, however, our goal is to obtain the flux 
on the circle, the canonical geometry for a finite absorber. 

2.4. Similarity Solutions for Semi-Infinite Leading Edges 

Before proceeding with our analysis, we mention a class of similarity solutions for "leading 
edges" which have relevance for the high-Pe limit of our problem. For this section only, we 
use the upper half plane, {haw > 0}, as our simple domain, fi^. For a straining velocity 
field, (piw) — Kew^, bringing fluid toward the plane, there is a classical similarity solution 
for the concentration profile, c{w,w) = erfc(-\/PeJ^Im w), for which the flux on the real 
axis is a constant, — 2-\/Pe|,/7r. (As usual, the complementary error function is defined 
byerfc(z) = (2/y^)/;°d^e-*^) 

For every conformal map, w = f{z), from the z plane to the upper half w plane, there 
corresponds another similarity solution ( [Cummings et aL\ 1999: [Bazant[ 2004), 

(j) = Ref{zf and c = erfc(v/Pe^Ini/(z)) for Im w > 0. (2.16) 

Note that the boundary condition V0 ^ x a,s \z\ ^ oo holds only if f{z) ~ ^/z as 
1 2; I 00, so these solutions correspond to more general flows near stagnation points. For 
the purposes of this paper, we discuss two choices for f{z): 

(a) f{z) = y/z, which maps the entire z plane, with the exception of the branch cut 
{Imz = 0,Re2: > 0}, onto Vtw The {x,y) coordinates of the z plane then correspond 
to the streamline coordinates, and the flux on the strip is (j{x) — ^jPe^Jirx. This well 
known formula can also be derived via replacing the upper limit of integration in 1)2. 9|l by 
00, and applying the Wiener-Hopf method of factorization UKreinl 1962; |Wijngaarden| 
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(a) (b) (c) 




Figure 3. The concentration profile (Pe = 1) calculated numerically in the interior of the 
circular disk is shown (a) in Cartesian coordinates, and (b) in polar coordinates. In (c) the 
numerically obtained values are shown for c(r, 6 = Q) (thin solid line) , c(r, 9 = tt) (thin dashed 
line), h{r,9 — 0) (thick solid line), and h{r,9 = n) (thick dashed line). 



1966: [Carrier et all 1983). This procedure is carried out systematically to all orders of 
approximation in Sec. 0] and via a different, rigorous method by [Margetis fc Choi| (2004). 

(6) f{z) = \fz + ^/\/z, which maps the exterior of the circular rim, {z : \z\ > 
1 and < argz < 2n}, onto IjBazantl 2004). This solution describes the advection- 
dominated (high-Pe) fixed point of the ADLA fractal-growth process (jBazant et al.l 
2003). From with z = re*^ and r > 1, we find 



c(r, 9) = erfc 



aie) 




(2.17) 



which is the leading-order solution for our BVP H2.12|l - H2.13|l as Fcq — Pe ^ oo. 

The fact that we consider finite absorbers leads to significant analytical and numerical 
complications, which are the focus of this paper. 



3. Numerical Solution 

3.1. Conformal Mapping to Polar Coordinates Inside the Unit Disk 

In this section we determine the concentration profile c by numerically solving the BVP 
H2.12|l - (|2.13l) . One of the difficulties in applying a numerical method is related to the fact 
that the region f2^ is unbounded. By invoking conformal invariance again, however, we 
can apply a transformation that leaves the BVP unchanged yet maps £7^, to a bounded 
region. In particular, we use inversion, g{w) = l/w, to map the disk exterior onto its 
interior. Physically, this corresponds to a dipole source of concentrated fiuid inside an 
absorbing circular cylinder, as shown in Fig. ISJa). With g{w) — re*^ the problem is 
expressed in the polar coordinates (r, 9) by 

+ {r^ + PcKl - r^) cos9}^ + r%^ + Pe(l + r^) sin^-^ = 0, (3.1a) 
or"' or 09'' o9 

c — 0, at r = 0, and c = 1, at r = 1, (3.1&) 
where r < 1 and < 6* < 2tt. A solution in the (r, 9) plane is shown in Fig. ^h). 

3.2. Analytical Treatment of Singularities 
Before we apply any numerical method to H3.1() directly, we note that the concentration 
profile c as a function of (r, 6*) exhibits singular behavior as r approaches 0. We first need 
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to modify H^.lafl in order to eliminate this behavior, which undermines the accuracy of our 
numerical method. From the similarity solution H2.17|l and Green's function of Sec. 12.21 
we obtain the leading-order behavior 



c(r, 9) — O < \/r exp 



Pe ( 2 r I sin 

r 



0. (3.2) 



First, the square-root limit c(r, 9 = 0) = 0{^/r) can hardly be dealt with in numerical 
methods because of the resulting diverging derivative near r = 0. Second, the essential 
singularity at r = 0, c(r, 9 = tt) = 0[^/r e^^'^/'') , forces c to change drastically near r = 0. 
Especially when Pe is large, this limiting behavior is extended even to r < 1 — 0(1/ \/Pe). 
To avoid this behavior, we define a function h(r, 9) by factoring out the leading-order 
singular behavior of c{r, 9) as 

^Jr exp jpe ( 2 - i - r ) sin^ ( ^ ) 1 h{r, 9) = r c{r, 9), (3.3) 



and apply the numerical method shown below directly to this h{r, 9). Note that the c is 
multiplied by r in the right-hand side of (|3.3|) to ensure that h{r, 9) = 0{r) as r ^ 0; 
thus, /i = at r = 0, which is the same condition as for c. Combining H3.1a|l and H3.3|l 
we obtain a PDF for h{r, 9): 

r'^+Pe{r-r')^+r^ + 2Per sm9^ + \Pe{r cos9 - 1) + h = 0. (3.4) 
or^ or o9'^ o9 I 4 J 

By comparison of (|3.1a|l and 1)3. 4|l . we note that the coefficients of the derivatives are 
simplified in (|3.4|l . Once h is determined from (|3.4|) . c is simply recovered via (|3.3|) . and 

a{9) is obtained as 



_ dh 
=i,s) 9r 



(r=l,9) ^ 

Fig. Ofc) shows how the singular behavior of c is mitigated by introducing the new 
variable h. 

3.3. Spectral Method 

The numerical differentiations with respect to the variables r and 9 are carried out by 
spectral methods (jTrefethenI 2000). The spatial nodes, the points {rj,9k) where the 
function is evaluated numerically, are determined by 

where j = 1, . . . , Nr and A; = 0, . . . , Ne; the nodes have higher density at the endpoints 
as shown in the frame of Fig.|3Jb). Once the function values are given at the nodes, the 
derivatives at these points are calculated by interpolation via Chebyshev's polynomials. 
This procedure is very efficient, as the error in the spectral method is known to decrease 
exponentially in the number of nodes UTrefethenl 2000). We used Nr = 50 and Ng = 100 
(practically Ng = 50 exploiting the symmetry in 9) for all the numerical results appearing 
in this paper. 

Fig. ^ shows the concentration profile, c{r,9), for Pe of different orders of magnitude. 
As Pe increases, there is an apparent crossover from a diffusion-dominated regime (a), 
where the concentration disturbance looks like a "cloud" extending in all directions, 
and an advection-dominated regime (c), where concentration gradients are confined to a 
narrow boundary layer which separates into a thin wake downstream. Understanding the 
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Figure 4. 



Concentration profiles for adsorption (or desorption) around tlie unit circular disk 
for (a) Pe = 0.01, (b) Pe = 1, and (c) Pe = 100. 
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for Pe = 0, 



The flux o"(^Q is plotted around the unit disk for different values of Pe: (a) a{6)/a{n) 
2~^, 2~*, 2~ , 1, 2^, CO, and (b) the values a{0 — tt) and a{9 = 0) for a range of Pe. 



crossover regime (b) to (c) is an important part of this paper, revisited below in Sec. [71 
following our analysis of the flux profile. 

Fig.lSlshows the flux density on the absorber, a{9), for different values of Pe. To check 
the validity of our numerical method, we compare the result for a with the asymptotic 
expansion for high Pe from Sec.Qlbelow, which can be numerically evaluated to any order. 
For this comparison we used the intermediate range of Peclet numbers 10~^ < Pe < 10^ 
for which the series converge fast enough yet the numerical method is stable. When 
Nr = 50 and Ng — 100 nodes are used for the discretization, the relative error measured 

as llcTnum - CTasymll / ||crasym|| is of the Order of 10~^ Or Smaller; here, CTnum and CTasym 

denote the numerical and the asymptotic solutions, respectively and the norm is defined 
as \ \a\\= max{6l e [0, 27r) : |cr(6l)|}. 

If Pe lies outside the given intermediate range, care should be exercised in using our 
numerical method as explained in the next subsection. For such cases, however, the 
asymptotic, analytical formulae of Sees. 0] and [S] give sufficently accurate solutions. The 
reader is referred to Sees. 01 and for details on the formulae and the comparisons with 
the numerical solution. 
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Figure 6. Contour plots of h for Pe — 100 (a) in the {r,8) plane, and (b) in the {f,0) plane. 



3.4. Adaptive Mesh for Very High Peclet Numbers 

A feature of the solution c that may undermine the accuracy of our numerics is the 
emergence of boundary layers for sufficiently large Pe. By virtue of ()4.1fcifl below we 
expect that for Pe ^ 1, 



dc 
dr 



ai0) 



(r=l, 



2^Pe/7r sin(6i/2) 
I/tt 



when 6'>0(l/\/Pc) 
when 6'<0(1/VPc). 



(3.5) 



So, c{r,9) has boundary layers near r = 1 and 9 ^ 0, 2tt whose widths are 0(l/vPe), as 
indicated in Fig. EJa). The layer aX 9 = corresponds to the "tail" shown in Fig.^Jc). 
Thus, the numerical method starts to break down when the node spacing becomes of the 
order of 1/ y/Pc. 

We next outline a technique to deal with the case of very high Pe within our numerical 
procedure. The idea is to introduce a set of independent variables, f = f(r) and 9 — 9{9), 
so that c is a sufficiently smooth function of f and 0. The following conditions are required 
for f(r) and 9{9): 

(3.6) 
(3.7) 

where the prime here denotes differentiation with respect to the argument. Once f{r) 
and 9{9) are defined, wc find the PDE for h{f,9) = h{r{r),9[9)) from l|3.4|l by applying 
the chain rule for the differentiations with respect to r and 9; for example. 



r(0) 0, f(l) = 1, f (1) O(VPe) 
9{Q) = 0, 9{2tt) = 27r, 9'{Q) = O(VPe), 



dh 
dr 



1 dh 
r'if) df ' 



d'^h 

Qj,2 



1 d^h r"(f) dh 



r'{fY dr^ r'{f)^ df 



We solve the resulting PDE numerically. A convenient choice for r(f) and 9(9) is 



r{f) 



sin(7rf), 9{9) = 9 



sin 9. 



(3.8) 



(3.9) 



The advantage of using f and 9 instead of r and 9 is illustrated in Fig.|Hlfor Pe — 100. The 
effects of the boundary layers in (r, 9) are notably suppressed in the formulation using 
(f, 9). In Secs.01and|Hlwe discuss the high- and low-Pe asymptotics and their comparisons 
with the solution determined numerically by the method of this section. 
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4. Direct Perturbation Analysis For "High" Peclet Numbers 

In this section, we derive an approximate analytical solution to the integral equation 
1)2.91) in terms of series expansions produced via suitable iterations in the coordinate 
space. We also obtain closed-form expressions for the terms of the iteration series as 
Pe-dependent multiple integrals. We show that the series is convergent for Pe > 0(1), 
and that retaining only a few of its terms produces accurate results even for Pe = 0(1). 
An iterative procedure in the Fourier domain that leads to the same results is given by 
[Margetis fc Choi| (2004) . 

4.1. Zeroth-order solution via the Wiener- Hopf method 

The starting point of the analysis is the observation that, as discussed in Sec. 12.41 the 
solution for the semi-infinite strip — 1 < a; < oo (in the variable notation of (|2.9|l ') 
provides the leading-order term of the high-Pe asymptotic expansion of the solution 
for the finite strip up to a distance 0(l/-\/Pe) from the endpoint a; = 1. In order to 
develop a systematical scheme for the correction terms, we symmetrize the kernel of 
()2.9(l and rescale the independent variable x using s = Pe{x + 1) while we define /x(s) 
by cr(x) = (v^/tt) Pee*/j,(s). (The factor {^/2/n)Pc is chosen for later convenience.) The 
integral equation H2.9|l thus becomes 

.2Pc 2 -s 

ds'Ko{\s- s'\)pl{s') ^ 0<s<2Pe. (4.1) 

v2 

An approximate solution fi ^ fiQ that is valid to the leading order in Pe is found by 
taking Pe — > oo in the upper limit of integration in 14.1|l . The resulting integral equation 
is 

ds'Ko{\s-s'\)tio{s')^^^^, 0<s. (4.2) 
V 2 

The solution /io(s) is obtained by the Wiener-Hopf technique UKreinl 1962: INoblel 

1988). Here we outline the basic steps of this method, which are also applied to other, 

similar integral equations below. First, we extend the validity of (|4.2|l to — oo < s < oo 

via modifying its right hand side, 

ds'Ko{\s - s'\) na{s') ^ —{e~''u{s) +p{s)}, -oo < s < oo, (4.3) 
V 2 

where fJ.ois) is taken to be zero for s < 0, u{s) is the Heaviside function (u(s) = for 
s < and u{s) = 1 for s > 0), and p{s) is an unknown function which has non-zero 
values only for s < 0. Next, we apply the Fourier transform in s to (|4.3|l . Defining the 
Fourier transform, /io(fc), of /Lto(s) as 

f°° f°° dk 

flo{k) = ds iio{s) e-''"' where ^o(s) = / — e'''' flo{k), (4.4) 



2tt 



yields 



7r/2o(fc) _ TT^ 



^ m 



1 + ik 



By simple algebraic manipulations the last equation becomes 



Mo(fc) 7I_ _ TT 

^/1 + ik 1 + ik^ ^ 



VI -ik-V2 
1 + ik 



+ Vl - ikp{k) 



(4.5) 



(4.6) 



where the left hand side defines a function analytic in the lower half k plane, Im k < e for 
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a small positive e, and the right hand side defines a function analytic in the upper half 
k plane, Imfc > — e; each of these functions vanishes as |fc| cx3 in the corresponding 
half plane. Thus, the two sides of H4.6|l together define an entire function of k, which is 
identically zero by Liouville's theorem IjCarrier et al.\ 1983). It follows that in the region 
of overlap, |Imfc| < e, the solution is /io(fc) = 7r(l + ik)~^/'^. Inversion of this formula 
yields 

/Xo(s) = y^^""' ^ <^ <°°- (4-7) 

4.2. Leading-order Uniformly Accurate Approximation 

The deviation of fJ-o{s) in (|4.7|l from the actual solution fi{s) of the finite strip is inter- 
preted as due to the effect of a fictitious, "misplaced" flux source lying in 2Pe < s, which 
is present in H4.2(l . Therefore a correction term must be found for /io(s) by placing a 
"correction source" on the original strip, < s < 2Pe, to compensate for the effect of the 
misplaced source. [Margetis fc Choi| (2004) further develop this approach and place it on 
a firm mathematical ground using Fourier transforms and a generalization of the Wiener- 
Hopf method. We proceed to calculate the correction to fiQ iteratively. Accordingly, the 
solution /x(s) is sought in terms of the series 

fl = fj.0 + fJ-i + fJ-2 + fJ-s + ■ ■■ + + • ■ ■ , (4-8) 

where each term, /i„(s), corresponds to suitable source corrections as described below. 

We consider a half-fine as the domain of the correction ^i(s), as we did to obtain 
but in the region —oo < s < 2Pe instead of the region < s < oo. The term ni is 
determined so that its effect compensates for the integrated effect of fJ.o{s) in the region 
2Pe < s. Hence, the correction /Lti(s) satisfies 

2Po (.oo 

ds'Ko{\s-s'\)^l,{s')^ ds'Ko{\s-s'\)fio{s'), s < 2Pe, (4.9) 

-oo J2Po 

where the right-hand side is known. The next-order corrections can be formulated and 
interpreted in a similar way; the correction fin compensates for the effect of the misplaced 
source corresponding to /i„_i, where fin-i and fin are defined in half- lines that together 
cover the entire real axis and overlap only in the region of the finite strip. In general, 
/x„>i satisfy the recursion relations 



Ko{\s-s'\)fin=2kis')ds' ^ / Ko{\s-s'\)fi^-i{s')ds', (4.10a) 

J —OO 

Ko{\v-v'\)fir,=2k+i{2Pe-v')dv' ^ / Ko{\v - v'\)fin-i{2Pe ~ v')dv\ (4.106) 



where we made the change of variable from s to v ~ 2Pe — s so that the integral equation 
for /i„=2fc+i has the same form as the one for fin=2k- The variables s and v are both 
positive (s > and v > 0), the left endpoint {x — —1) of the strip corresponds to s = 
and the right endpoint (x = 1) corresponds to w = 0. 

The various /i„ can be obtained successively, order by order, by applying the Wiener- 
Hopf method IjKreinl 1962: INobleL 1988) directly to H4.10|) . but the procedure becomes 
increasingly cumbersome with n. Instead, we propose a systematic procedure that fa- 
cilitates the derivation of a closed- form expression for each fin- For this purpose, we 
introduce an operator, £, that relates fin-i and fin by £[/i„_i] = fin- By H4.1()|l C is 
linear. In order to obtain fii, we notice that the leading-order solution fio{s) can be 
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represented as an integral over a variable, to, 

A^o^^/^e-^^ r dioe~^(i+*o). (4.II) 



s 

Then C acts on /ig to yield /ii as f 



^ dto C[e-<'+'<>^], (4.12) 



where the order of C and integration is safely interchanged. The advantage of using the 
to-representation is that e^*^^~'"*o-', as a function of s, has a Fourier transform simpler 
than the Fourier transform of /io('S) itself. The function £[e~^'^^+*o)] is found by the 
Wiener-Hopf method (Krein, 1962) as described in Sec. 14.21 and /ii follows by 14.12|l : 



e-2Pe(l+to) 



^ e-'" - Tr^2 + tl e''(i+*o)erfcy't.(2 + t2) 



(4.13) 



Aii(w) = ifo(2Pe)-^ - / dto e-(2Pc-«)(i+t^)erfc./i;(2 + t2). (4.14) 
^/7r?; / V 



Because and /ii are accurate arbitrarily close to the left edge (s = 0) and the right 
edge {v = 0) of the strip, respectively, /xq + /ii yields a leading-order approximation for 
;U as Pe ^ 00 valid over the entire finite strip. The corresponding approximation for the 
flux (T = ('\/2/7r)Pce*/i(s) is cr ^ ui + (T2 , which is given by 



,(hi),^^_o./Pe J 1 , Xo(2Pe)e2P- 



TT [ ^2(1 + 2;) t:^2{1-x) 



- dT ^^erfcV(2Pc + t^){1 - 

J \/2¥ 

for the geometry of the finite strip, and 



(4.15) 

X) \ 



a{9) - (7('")(0) = 2a/— ( I sin -I + -ifo(2Pe)e' 

V TT 1^ 2 TT 



2Pccose. Qg_| 



/27r 



J dT e-(i+cose)r^g^f^^(2Pe + t^){1 - cost 



(4.16) 



for the geometry of the unit circular disk. In the last formula we changed the variable to 
r — \/Pe to. An expansion similar to H4.15|l has been obtained by |Chugunov fc Kornev| 
(1986) and |Kornev fc Chugunov| (1988) in the physical context of artificial freezing |; 
these authors did not use the operator C in their method and apparently did not obtain 
higher-order terms. An elaborate mathematical procedure for asymptotic solutions to 
the relevant class of integral equations on the basis of kernel approximations is described 
in IjAleksandrov & BelokonI 1968) and l|Aleksandrov fc PozharskUl 1999). We note that 
ao{x) and ai(x) are singular at x = — 1 and x = 1, respectively, the edges of the finite 
strip; these singularities are properly removed after the strip is mapped onto the unit 
disk. 

f In the remaining integrals of this section the range of integration is understood to be from 
—00 to 00 unless it is stated other wise. 

J We could not verify whet her (14.1511 is equivalent to (13) in [Chugunov fc Kornev| (1986) or 
(3.7) in |Kornev fc Chugunov| (1988). 
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4.3. Exact Higher-Order Terms 

In an effort to obtain further insight into the nature of the solution a of (|2.9I) , we next 
derive exact, closed-form expressions for ^„ for all n in terms of iterated, multiple in- 
tegrals. For this purpose, we exploit the C operator introduced above. We observe that 
H4.13|l has the integral representation 



dti' 



-2Pc(l+t^) 



7r./2T4 2 + t', + t 



(4.17) 



where the last factor in the integrand has the same form as the e~*'-^"'"'o) term, on which 
the C acts in H4.12|l . with the s being replaced by v and the to being replaced by ti. It 
follows that the Hn is expressed as the n-th power of £ acting on a term of the form 
g-s(i+t )^ £njg-s(i+t Thus, by induction, ^„ is expressed as an iterated, multiple 
integral of the independent variable s or v, 



-2nPc 



dto dtiiQoRi) dt2{QiR2)-- - /dt„(g„_ii?„)e-"(i+*"\ (4.18) 



(4.19) 



where u=s for even n and u = v for odd n, and Qn and i?„ are defined by 



-2Pct„ 



Qr. 



R„ 



t 



2 + <Li 



tl 



^V2 + tl 

As indicated from 14.13|) . /x„(u) has has the singularity at u = which comes 

from the last integral of (|4.18() . Thus we can single out the singular behavior as fin (u) = 
y^TrJu e~2"P°~" Fn{u) and the singular-free part, Fn, is given by 



dtr 



Fn{u) = 

■J dtn-l{Qn-2Rn-l) 



dtiiQoRi) J dt2{QiR2) ■ ■ ■ 

Qn-i - e2-'(^P<=-"'*"-erfcy^«(2 + tl_,) 



(4.20) 



where (|4.13|l and 1)4. 17|) are used to evaluate the last integral. 
Thus, the nth term for the flux on the strip, cr„(a;), is given by 



Pe 



-2nPo 



crn=2kix) = 2y — 

0'n=2fe+l(a;) = 



^2(1 +X) 

Pe g-2(n-x)Pc 



F„(Pe(l + a;)) 
i^„(Pe(l-a;)). 



(4.21a) 
(4.21&) 

(4.22a) 
(4.226) 



TT ^2{1+X) 

On the unit circle cr„((?) is free of singularities in 6 and is expressed as 
cT„=2fc(^^) = 2J— e-2"P°|sin-| F„(Pe(l-t-cos^)), 

V TT Z 

a„=2fc+i(0) = 2^/^e-2(»-cos9)Pe| ^1 F„(Pe(l - cos9) ) . 

V TT Z 

It has not been possible to evaluate (t„ from H4.22|l in simple closed form, except for 
n = 1, 2, as described by H4.15|l . However, the numerical integrations over the variables 
tj{j = \,2,...n) for each ct„ can be carried out efficiently by using recursion. 

We verify that the sum for the ffux, a <tq + ui + . . . -\- am calculated for finite n via 
the numerical integration of 14.22|l . indeed approaches the numerical solution of Sec.|3l 
In Fig. 13 we show a comparison of the numerical solution of Sec. 01 for (7(0) and a{'n) 
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(a) (b) 




Figure 7. Asymptotic approximations for high Pe compared with our numerical solution: 
(a) Upstream flux, a(9 = 7r;Pe) and (b) downstream flux, (t{9 — 0;Pe). The integrals in the 
asymptotic corrections, (t„ (1 < n < 5), were performed numerically. 



with formula H4.22|l for different values of n. Because <Jn=2k(P) and (Jn=2k+i(Q) vanish 
&% 9 = Q and 6 — n, respectively, the term an=2k{0) affects only 17(77) whereas the term 
o'n=2fe+i(^) affects only ct(0). 

Remarkably, with only a few terms, our approximation is uniformly accurate down to 
values of Pe of order 10~^ or lower, which could hardly be called "high", while some 
correction terms /z„>2 may not be small. This behavior suggests that there may be 
an intermediate region of overlap between asymptotic approximations for high and low 
Pe. Indeed, by combining such approximations below, we will construct a very accurate 
approximation for all 9 and all Pe. 

The closed-form expression of (T„ also serves as another "numerical method" for high 
Pe. The multiple integrals in (14.20(1 can be numerically evaluated in a recursive way; 
the intermediate calculation steps for Fk recur in the calculation for all F^^k- Thus, the 
computational cost for X]fc=o '^^ same as that for (t„, which scales linearly with n. 

4.4. Convergence of the Asymptotic Series 

We next discuss the convergence properties of the asymptotic series ^'^^ 

the unit circle. Because the maximum of ct„(0) occurs at = tt for n = 2k and ai 9 = 

for ri = 2fc + 1, it is readily shown that 

||,.||.275 and I ll--*ll-;:r,f;»"IW,, ,4,23) 

From (|4.2U|I . ||F„|| is simply — 0). In particular, for n = 1, 2, ||i^„|| are evaluated 

in simpler forms: 

||Fi|| = ^i^o(2Pe), \\F2\\^'-^{^^^^^^- TdtKoitA. (4.24) 



TT TT^ 



2Pe 



We have not been able to further simplify the expressions for ||-F'n>2||- We now show 
that each ||F„|| is bounded by a function of Pe that ensures convergence of the series 
J2n <^n for Pe > 0(1). By noting that g„ < e-^P^*' /V2 and i?„ < I/tt, we find 

IIKII <7r-"/2(4Pe)-"/2, (4.25) 

for any Pe > 0. Hence, by H4.23|l the series is characterized by geometric con- 

vergence in the parameter Pe[~^/^] for Pe > 0(1). Finally, for large Pe and fixed n 
the asymptotic behavior of \\Fn\\ with respect to Pe is obtained via scaling the original 



(a) 



Advection-dijfusion in two dimensions 
(b) 



17 





A 


A 


A 


A' 


A 


A 


.A' 


.A.- 


A 


A 


A 


A 































□ 
























A 














o 















□ 



































o 


o 







■0 









o 


« 




o 















A 


Pe = 


1 






O 
















□ 


Pe = 


10"' 


























Pe = 


10"- 

























Pc = 


10"' 









2 




4 




6 




8 




10 




12 



10° 4i-^:;;A--A...A...A...A .a...a.. a..,a...a.,a. 
5 ° □ 

«... , . 



^ Pe=10' 
° Pe=10" 
Pe= 1 



'A...A...A...J1 



8 10 

n 




Figure 8. Numerical evidence for the convergence of the iteration series o"„ ~ a. (a) 
The ratio ||F„ 1 1/| |F„+i 1 1 as a function of n for different values of Pe. The convergence of 
X^n guaranteed if | |/| 1 1 > 1 by virtue of l|4.23^ . (b) The term ||f„|| as a func- 
tion of n for different values of Pe. (c) The decay rate p(Pe) of \\F„\\ as a function of Pe, 
where p = lim„_,oo(| l^n 1 1/| |Jn+i ||)- The asymptotic behavior H4.26|l is shown to be attained for 
Pe > 10. (d) The ratio of X]fc=i '^k{d) to (t(6) as a function of small values of Pe and 6^ = 0, tt, 
where (j{6) is evaluated numerically by the method of Sec.|3 

,n- 1), 



variables tk as — tk v 2Pe (fc = 0, . 

MM 1 / 

Fn>l" ' 



1 



as Pe oo. (4.26) 



Vl6\/^PeP/2], 

Formulae H4.25|l and H4.26(l indicate that the series ^" converges geometrically for a 
wide range of Pe. 

We check numerically that ||i^n|| decays exponentially in n for a wide range of Pe, 
\\Pn\\ ~ as n ^ OO where p = p{Pe) > is the "decay rate" of \\Fn\\, which is 
independent of n. For this purpose, we examine the ratio ||F„||/||i^„+i|| as a function of 
both n and Pe, expecting that this ratio approaches a constant for fixed Pe as n becomes 
sufficiently large, as shown in Figs.|S^a), (b), (c). From Fig.lSJd) we find that the relative 
error in the approximation of a by the sum X]fe=i "'fe becomes higher than 1% only when 
Pe < 6.5 X 10-3. 



5. Uniformly Accurate Asymptotics For Low Peclet Numbers 

In this section we solve approximately the integral equation H2.9|l for the surface fiux 
(j{x), for all X in (—1, 1), when Pe is sufficiently small, Pe < 0(1). For this purpose, we 
define the dependent variable if{x) — 6^^°^ a{x). Equation (|2.9I) thus becomes 

J dx' Ko{Pe\x-x'\)(p{x') = Tre-'^'''', -1< a; < 1. (5.1) 
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Because the argument of the kernel is also sufficiently small, Pe|x — a;'| < 0(1), we invoke 
the expansion 



K^{Ve\x-x'\) ^ -/o(Pe|x-x'|) \J^\x-x'\\ + 



(m!)2 

m=0 ^ ' 



, (5.2) 



where it is understood that 



/o(Pe|x - x'l) ^ 7zk^ (5-3) 

m=0 ^ 

and ^^{z) is the logarithmic derivative of the Gamma function, 'i/;(z) = ^lnr(z). It 
was first pointed out by IPearsonI (1957) that the resulting integral equation can be 
solved exactly for any finite number of terms, M , but the procedure becomes increasingly 
cumbersome with M . 

To leading order we consider M = in H5.2|l and H5.3|l . Equation H5.1|) thus reduces to 
a variant of Carleman's equation l|Carlemanl 1922), 

dx' ln(|a; - x'\) ^q{x') = Ci ~ ttc'^"'' , (5.4) 



where 7 = — ■0(1) — 0.577215 • • • is the Euler number, ^pq{x) ^ ^{x) is the corresponding 
approximation for tf{x), and Ci is the constant 

Ci = -[7 + ln(Pe/2)] / dx' ipo{x). (5.5) 



J -I 

Following [Carrier et aLl Cl983'). we introduce the complex function 



$(,) ^ ^^ZH /' dx' ^ (5.6) 



2m 

and single out the limit values 

/i 2 /"iT^^ ( f\ 

<^^(x) = \im^x±ie)=± ^ lim / dz' (5.7) 

e^o 27r e^o7_i^,^ z'-x 

by which the integral equation (j5.4|l is equivalent to the equations: 



d 



<i>+{x) - ^-{x) ^ 4- / d^' H\x - ^'\) M^') = -PtVl - e-P°^ (5.8) 

TT dx 



(5.9) 



$+(a;) + $~(a;) = Vl - Resi '^"'-^ ^ ; z' = xl = i ^1 - j;2 ipo(x). 

{z' — X J 

First, we find ^{z) via applying directly the Mittag-Leffler expansion theorem to H5.8|l 
HCarrier et aU 1983): 



X 



$(z) = --^ I Pe / dx' " ^, 6-^^=="' + A , (5.10) 
27ri \ a; — z / 

where — (27ri)^^ A is the limit as 2 ^ 00 of the function $(2:); by inspection of (|5.6|l . 

^ = y da;Vo(a;')- (5.11) 
We recognize that the constant Ci of H5.5|l is Ci = — [7 + ln(Pe/2)]yl. 
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Next, we obtain the approximate solution (po{x) in terms of this A by 1)5. 8|l : 



ipo{x) 



1 



TT\/l — 



Pe [P) / dx 



a; — a; 
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(5.12) 



where (P) J^j^ denotes the principal value of the integral. In order to determine the 
unknown constant A, we multiply both sides of (|5.4|l by (1 — a;^)^^/^ and integrate over 
(—1, 1) by use of the elementary integral ( Carrier et al.i 1983) 



^^^ln(|.-.'|) 



-7rln2. 



(5.13) 



A few comments on this result are in order. It can be obtained via differentiating the left- 
hand side with respect to a;', and thus converting the integral to a Cauchy principal value 
which is found directly to vanish identically. Hence, the original integral is independent 
of x' and can be evaluated for a;' = by changing the variable to a; = (^ — l/^)/(2«), 
where ^ moves on the unit circle, and applying the residue theorem IjCarrier "eFoIl 1983). 
We thus find 
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It is straightforward to carry out the integrations in H5.15II . The second integral on the 
right-hand side is simply a modified Bessel function of the first kind: 

"1 ^-Pcx 



dx 



dtt 



-Pc COS t 



= 7rJo(iPe) = 7r/o(Pe). 



(5.16) 



The remaining integral can be converted to one that is directly amenable to numerical 
computation for Pe < 0(1). By defining 



I(Pe;x) = (P) dx' ^ ' e-P°("'-") 



a;' — a; 



we evaluate the derivative 
dl 



-Pea; 



a;'2 e-P°^' 
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where Ii, is the modified Bessel function of the first kind. An expression for the integral 
H5.17|l then follows by direct integration in Pc of (|5.18|l : 



I(Pe; x) = X(0; x) - n f dt e*^ 

Jo i 



where 



Hence, 



I(Q;x) = -(1 - a;2) 



I(Pe; x) 



dx 



, ln(|a;-a;^|) 

dx — , Tlx 

VI - x'-^ 



Po 



dt e 



(5.19) 
(5.20) 
(5.21) 



20 

(a) 



J. Choi, D. Margetis, T. M. Squires, and M. Z. Bazant 

(b) 



3 
2.5 
2 

1 

0.5 



Pe 



Pe 



Figure 9. The asymptotic approximation ct^'"-' (0; Pe) in l|5.24^ compared to the numerical 
solution of Sec.Elfor a range of low Pe. (a) Upstream flux, a^'^°\9 = tt; Pe), and (b) downstream 
flux, cr''°^(e = 0;Pe). 



The approximation 

7 + ln(Pe/4) ~ -/fo(Pe/2), (5.22) 
which becomes useful in Sec. where we construct a uniform formula for all Pe and local 



coordinate of the boundary, and the use of H5.16(l . H5.21(l and a{z 
low-Pe approximation for the flux on the boundary of the finite strip, a 
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Hence, by virtue of H2.15|l . the flux on the unit circle is 
,(lo)..^_ ^o(Pc) 



HO) 



.,Poc 



A'o(Pc/2) 



Pe 



Pc 



dt e* 



hit) 



(5.24) 



We note in passing that the integral in (|5 
Pe. With the series expansions 



can be expressed as a powers series in 
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it is straightforward to derive the expansion 
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A few comments on formula ()5.23|) are in order. lAleksandrov &: 13eiokonl (1967) ex- 
panded to high orders the kernels of the relevant class of singular integral equations and 
derived in generality a more accurate yet complicated formula for the solution. The pro- 
cedure here, though being based on simply taking M = in (|5.2|) and H5.3|l , leaves intact 
the right-hand side of H5.1|l and applies (|5.22() : our approximate formula for aix) turns 
out to be accurate for an extended range of low Pe. Fig. shows a comparison of H5.24|l 
with the numerical solution of Sec.|31for a range of low Peclet numbers. 
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Figure 10. Plots of formulae ijCT))! for a^*"'' and if^^^ for cr^'"' versus Pe in the same graph 
as the plot for the solution o evaluated numerically by the method of Sec^ (a) Upstream flux, 
a{9 = 0; Pe). (b) Downstream flux, o{9 — tv; Pe). 



6. Uniformly Accurate Formula for All Positions and Peclet Numbers 

6.1. Connecting the High and Low Pe Approximations for the flux 

In Sees. 01 and 13 we derived asymptotic formulae for the surface flux a on the boundary 
of the unit circular disk (or finite strip) that are valid for sufficiently high or sufficiently 
low Pe; these expressions, ct*^^') and cr('°' in H4.16(l and l|5.24|l . respectively, hold for all 
values of the local coordinate of the absorber boundary, although we did not analyze to 
what extent the approximations are uniformly valid in a rigorous mathematical sense. 
Comparisons with the numerical results in Figs. [7| and |^ show that the approximations 
are comparably accurate for the rear stagnation point [9 — 0) and forward stagnation 
point {9 = tt). We have also checked that, for fixed Pe, the errors are also comparable 
at intermediate values of the local coordinate {0 < 9 < n). In Fig. EH we show that the 
two approximations, ct^'"-* and ct^'°^ , nearly overlap for some range of values Pe near the 
value Pe = 0.1, while remaining remarkably close to the "exact" numerical solution from 
Sec.El 

The existence of a regime of overlapping accuracy allows us to construct an analytical 
formula for a, accurate for all values of Pe and the local coordinate of the absorbing 
boundary, 9, by smoothly connecting cr*^^'^ and ct*^'"-*. A similar overlapping accuracy 
was found for a three-dimensional problem of heat or mass transfer in a steady shear 
flow by lPhillips] (1990), who combined high- and low-Pe expansions only for the Nusselt 
number, Nu, using singular perturbation. The dependence of the flux a on Pe can be 
described heuristically by a formula of the form 

a{9;Pe) - ct^™"") = a^^'\9;Pe)U(Pe/Po) + (T'^^°\9;Pe) [1 - U{Pe/Po)], (6.1) 

for the entire range of Pe and 9; U{x) is a family of smooth functions defined for x = 
Pe/Pp > that at least satisfy the conditions 

lim U{x) = 0, lim U{x) = 1- (6.2) 

A simple choice for the step-like function, which yields a rather accurate formula for 
cr, is 

Zi(x;a) = e^/(i-""°\ a>0. (6.3) 

We note that there are two free parameters in 1)6. 3|l with x = Pg/Pq: a and Pq. The 
parameter Pq corresponds to a value of Pe in the region of overlap of formulae H4.16|l 
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Figure 11. Plots of the connection formula, o-''^™"^ from Eqs. I|6.1^ and \iS.'A\ . in comparison 
to the numerical solution of Sec. |3 (a) Upstream flux, a[9 = 7r;Pe), and (b) downstream flux, 
(j[0 = tt; Pe). The relative error of a'™""' compared to the numerical solution is less than 1.75% 
for aU values of Pe and 6, and it becomes negligibly small at high and low Pe for all 6. 



and H5.24(l : in principle, Pq may depend on the local coordinate, which is 6 for the unit 
circle. The parameter a determines the steepness of the curve hi{x) near x = 1- 

We find that a good fit with the numerical solution is achieved for a = 2 and Pq = 1/6, 
as shown in Fig. 1111 where (7*^'^°""^(0) is plotted versus Pe for 6* = tt (upstream flux) and 
9 — (downstream flux). The relative error is less than 1.75% for all values of Pe and 6. 
Because we have the exact Green's function (|2.7(l . this uniform accuracy also carries over 
to the solution of the entire concentration field, c(a;,y;Pe), obtained from the integral 



Our analytical approximation also describes an absorber of arbitrary shape obtained 
by conformal mapping, z = g{w) {w = e*^). The fiux distribution on the unit circle, 
p.(conn) ^^-j ^ is transformed to the new surface using (|2.14|l . Because the fiux is propor- 
tional to a gradient, it is locally amplified by a factor of |(7'(ti;)|^^, which may cause 
relative errors larger than 1.75% in some locations, e.g. near a cusp, where conformality 
breaks down (g' = 0). For a well behaved univalent mapping, however, the approximation 
should remain very accurate for all positions, ^(e*^) and Pe, so the general BVP may be 
considered effectively solved. 

6.2. The Total Flux to the Absorber 

It is straightforward to obtain a uniformly accurate approximation of the Nusselt number 
by integrating the fiux on the unit circle or the finite strip: 

/'27T /• 1 

Nu(Pc)=/ ddcr{9;Pc) ^2 dxa{x;Pe). (6.4) 
Jo J~i 

The leading-order expressions of Nu(Pe) in the high and low Pe limits are thus obtained by 
integrating (j4.15|l and l|5.23|l using integration by parts and the identity xIq (x) — [xli (x)]' 
HGradshteyn fc Ryzhik| 1980), 

Nu(''')(Pe) = ^/^e-2P^i^o(2Pe) + Pee^P^erf (2\/Pc)[ifo(2Pc) + i^i(2Pe)] } , (6.5a) 



Nu('°)(Pe) = 27: {;^^|^ + (Pe)'[/i(Pe)' - /o(Pe)^] + Pe/o(Pe) /i(Pe)| . (6.5&) 
The uniform analytical approximation for Nu(Pe) follows as 

Nu(™-)(Pe) = Nu('^')(Pe) e^/^i-^""''') + Nu('°)(Pe) [1 - e^/^i-^""''')]. (6.6) 
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Figure 12. The Nusselt number, Nu, which gives the total flux to the absorber, versus the Peclet 
number, quantilying the importance of advection compared to diffusion. In (a) the asymptotic 
expressions l|6.5^ overlap for 0.1 < Pe < 0.4. In (b) the analytical connection formula l|6.6^ 
compares very well with the "exact" numerical result from Sec.|3 The results hold for absorbers 
of arbitrary shape, as long as Pe is the renormalized Peclet number, Pe = A\ Peo, where A\ is 
the conformal radius and Peo is the bare Peclet number for the unit circle. 

As shown in Fig. 1121 this analytical result is quite accurate over the entire range of Peclet 
numbers, and it becomes exact as Pe — > and Pe — s- oo. The maximum relative error 
is found to be 0.53%. We are not aware of any such analytical formula for the complete 
Nusselt-Peclet relation of a finite absorber. 

It may come as a surprise that the same result (|6.5ll also holds for an absorber of 
arbitrary shape, obtained by conformal mapping of the unit circle, z — g{'w), without 
any additional error. The reason is that the total flux through any contour is preserved 
exactly under every conformal mapping of a conformally invariant BVP (|Bazant| 2004). 
When computing Nu from (|6.4() for another shape, therefore, one must simply use the 
renormalized Peclet number, Pe = ^iPco, equal to the conformal radius, Ai, times the 
bare Peclet number, Pco, for the unit circle. 



7. Conclusion 

7.1. Summary of Results 

We have performed a detailed study of the BVP (|1.3|l - (|1.5(l for a finite absorber of arbi- 
trary cross section in a steady two-dimensional potential flow. Our focus has been the flux 
profile on the boundary, u, from which the concentration can be obtain everywhere in 
the plane by convolution with the known Green's function. We have explicitly considered 
several simple cases, notably the canonical problem of uniform flow past a circular cylin- 
der, which can be mapped to arbitrary geometries by conformal mapping, as described 
in Sec. [2 

In Sec. 13 we presented an efficient numerical method to solve the BVP. We began by 
mapping the BVP to the inside of the unit circle in order to work with a bounded domain. 
We then eliminated some "far-field" singularities (at the origin) using exact asymptotics 
and applied a spectral method in polar coordinates, with exponential convergence in the 
number of nodes. We also used an adaptive mesh to deal with boundary layers at very 
high Peclet numbers. The results, which are illustrated in Figs.|2]and0I provided reliable 
tests of our analytical approximations. 

In Sees. 0] and m we derived distinct asymptotic expansions for the flux on the bound- 
ary of the finite strip for high and low Pe, respectively, starting from a well known 
integral equation in streamline coordinates ( |Wijngaarden| 1966), which has been studied 



24 



J. Choi, D. Margetis, T. M. Squires, and M. Z. Bazant 



extensively in the theory of sohdification and freezing in a flowing melt fMaksimov' 1976; 
Chugunov fc KornevJ 198 6;[Kome v fc Chu gunovj 1988;IKornev fc MukhamaduUina. 1994; 
XTmio^t g/.l 1994: lAlimov et al.\ 1998; [Cummings fc Kornev[ 1999; [Cummings et al.\ 
1999). We used some original variations on classical techniques from the theory of sin- 
gular integral equations, including the Wiener-Hopf method of factorization, to obtain 
improvements on previous approximations. 

Noteworthy features of our expansions for high Pe are: (i) The summands admit closed- 
form expressions in terms of multiple integrals l|4.22|l . which are straightforward to eval- 
uate numerically; (ii) the expansions converge for all distances along the strip for a wide 
range of Pe, Pe > O(10~^); and (iii) only the small number of terms in H4.16|l can be 
retained for reasonable accuracy though the number increases as Pe decreases. On the 
other hand, our asymptotic formula H5.24|l for low Pe is accurate for Pe < 10"^, which 
renders it possible to have a region where the two expansions overlap. 

Therefore, we were able to construct an ad hoc analytical connection formula (|6.1|l . 
which reproduces our numerical results for the flux to a circular absorber with less than 
1.75% relative error for all angles and Peclet numbers, as shown in Fig. ^2 We also 
predicted the Nu(Pe) relation Ht).4|l with comparable accuracy, as shown in Fig.^l These 
results constitute a nearly exact, analytical solution to the BVP (|1.3|I - H1.5|I . 

The main contribution of our work is a unified description of the crossover regime, 
intepolating between the well known asymptotic limits of high and low Pe. As such, 
we can draw some mathematical and physical conclusions about the transition in the 
following sections. 



7.2. Mathematical Discussion 

A posteriori, we may try to understand why perturbation methods, which ostensibly 
require extreme values of Pe, produce a very accurate analytical solution for all values of 
Pe. One technical reason is that our high-Pe approximation, cr^*^'), is not the usual asymp- 
totic series of singular perturbation, since we have essentially "summed" parts of a naive 
expansion exactly in the Bessel function terms of (|4.16() . This allows the approximation 
to be more accurate than a simple powcr-scrics type of expansion, presumably extending 
its validity to somewhat higher Pe. Still, the neglected higher-order terms involve powers 
of Pe, which become important as Pe becomes small. 

As we have mentioned throughout the paper, similarity solutions govern the asymptotic 
limits. For Pe 0, we are perturbing around the trivial similarity solution with uniform 
diffusive flux from infinity, a =constant and c{r,9) oc Inr. As shown in Fig ^a) for 
Pe — 0.01, a small amount of fluid flow changes this picture only slightly near the disk, 
by favoring flux to one side compared to the other. When Pe = 1, as in Fig. Efb), the 
region of depleted concentration begins to be swept significantly downstream by the flow. 
This causes the low-Pe approximation to break down near the rear stagnation point, while 
remaining fairly accurate near the forward stagnation point, as shown in Fig.|3 

The high-Pe approximation is derived by perturbing around a different similarity so- 
lution (|2.17|l for an absorbing circular rim on an absorbing flat plate HBazantl 2004). 
In this advection-dominated regime, there is a thin diffusion layer of width, 0(l/-\/Pe), 
around the disk, as shown in Fig^fc), which provides the first term in the approximation 
l|5.24|l . The asymptotic corrections in Sec. 01 are obtained by effectively removing the 
"false plate" from the similarity solution, downstream from the disk. 

To understand the influence of downstream perturbations, consider the Green's tunc- 
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tion (|2.7|l . which has the asymptotic form 

gPo(K-r) 

G{x,y) ^ —j=== as r ^ y + y"^ ~* oo. (7-1) 

The Green's function decays exponentially at the scale of Pe in all directions, except 
precisely downstream, where it is long-rang ed: G{x,Q) - l/V87rPea; as a; — > oo for y — 0. 
Therefore, all corrections to the leading-order similarity solution 12.17|l decay exponen- 
tially upstream beyond an 0(Pe) distance from the rear stagnation point. Our high-Pe 
approximation, ctC"), in (|4.1()l) captures the first such correction, which is needed for uni- 
form accuracy over the disk. Further corrections are 0{e~'^^°), as is clear from formulae 

The fact that the approximation breaks down when e^^°, rather than Pe, gets close 
to unity explains the fortuitous accuracy of cr*^'''-' down to Pe = 0.1. Higher-order terms 
further extend the region of accuracy by orders of magnitude, e.g., to Pe — 0(10^'^) for 
five terms, as shown in Fig.[71 Because the approximation is valid for a wide range of Pe 
that would not be considered "high" , it overlaps with the low-Pe approximation. 

We also mention some directions for further analysis. The analytical treatment on the 
basis of an integral equation essentially avoids the complications of boundary-layer theory 
applied directly to the BVP fHinch' 1991). In that sense, it resembles renormalization- 
group (RG) methods for PDEs ( Goldcnfcld , 1992: [CIienl 1996), which provide a general 
means of deriving asymptotic expansions in place of traditional problem-specific singular- 
perturbation methods. An attractive feature of RG methods is that they promise to 
produce globally valid approximations directly, without the need to combine distinct 
overlapping expansions, as we have done. It would be interesting to see if this could 
possibly be accomplished for our BVP, and, if so, whether the resulting approximations 
are any simpler or more accurate than ours. We leave this as an open challenge to RG 
aficionados. 

From the point of view of mathematical methods, another interesting observation is 
that the BVP (|1.3|I - H1.5|I can be formulated in analogous way to the problem of wave 
scattering by a finite strip ( |Myers[ 1965) in acoustics or electromagnetics, which has 
some variants with other geometries such as the scattering by a broken corner ( |Myers[ 
1984). Such problems, where the requisite PDEs are or can be reduced to Helmholtz- 
type equations by simple transformations, can be described alternatively by finite sets of 
nonlinear ODEs in which the length of the strip, or the Peclet number Pe in the present 
case, is the independent variable. The underlying method is an improvement over Latta's 
method l|Latta[ 1956) for the solution of a class of singular integral equations via their 
exact conversion to ODEs. 

Another technical question is how to place the method of high-Pe expansion pursued 
in Sec. 0] on a more firm mathematical ground. To address this issue, two of us have 
developed an equivalent approach based on a generalization of the Wiener-Hopf technique 
applied to the BVP (|1.3|l - l|1.5|l in the Fourier domain, which yields the same results as 
the iteration scheme of Sec. ^ | |Margetis fc Choi| 2004). Since the Wiener-Hopf method 
is well known in this context for a semi- infinite strip ({Carrier et aZ.I 1983), it would be 
interesting to ask how the method might be extended for the present case of a finite strip, 
or perhaps more general situations, such as multiple strips, or three-dimensional objects. 

A more careful comparison should be made of our numerical solutions with the works by 
QCornev fc Mukham adullina., (1994). Our method in Sec. 13 seems to require more compu- 
tational resources to obtain the same level of accuracy because it solves a two-dimensional 
BVP, whereas the other method solves a one-dimensional integral equation. Aspects of 
the technique in Sec. 14.31 also need to be studied further, especially the number of the 
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Figure 13. The locus of points, for various values of Pe, (solid lines) where the concentration c 
attains its maximum along the streamlines (dashed lines) around a circular disk for desorption 
into an unconcentrated fluid. (In the equivalent problem of absorption from a concentrated fluid, 
the solid curves give the minimum concentration along streamlines for different Pe. ) 

required terms versus Pe, the accuracy for Pe ^ 1 and the propagation of error through 
the requisite multiple integrals. 

7.3. Physical Discussion 

Our analysis provides quantitative information about the crossover between diffusion- 
dominated and advection-doniinated adsorption, which is important in applications. As 
described above, the former regime, at low Pe, is characterized by a broad diffusion layer 
extending in all directions like a "cloud", with the flow causing only a minor broken 
symmetry, as shown in Fig. ^a). In contrast, at high Pe, concentration gradients are 
confined to a thin boundary layer around the object, which separates into a long thin 
"wake" near the rear stagnation point, as shown in Fig.0Jc). 

What is the critical Pe for the cloud-to- wake transition? Of course, the answer depends 
on one's definition, but the shape of the object also plays a role. In streamline coordinates, 
the transition is not so obvious since the object corresponds to an infinitely thin strip, 
which cannot "shield" a thinner, finite-sized wake. Objects of finite thickness, however, 
do show a clear transition. A uniformly valid formula for c{x, y, Pe) in any geometry may 
be obtained by inserting our expression for the flux to the strip into the convolution 
integral (|2.8|) and conformally mapping from the desired shape to the strip, w — f{z). 

The cloud-to- wake transition is related to f{z), as we illustrate for the case of the 
disk, f{z) = z + 1/z. The analytical approximation is nearly indistinguishable from 
the numerical solution in Fig. ^ so we consider the latter. In Fig. 1131 we show curves 
where, for different values of Pe, the concentration change, relative to the background, is 
maximal along streamlines. Far from the disk, they approach parabolae, Pey^ = Ax, due 
to balance between diffusion in the y-direction scaling as y ~ and advection in the 
a;-direction scaling as a; ~ Vet. Near the object, the map causes a significant distortion 
of the curves, which allows us to define the critical value, Pe = 60, above which they 
become non-monotonic functions of x. At larger values of Pe, the curves get sucked back 
into a thinner wake region behind the disk, signifying the dominance of advection. 

It is interesting to note the limiting wake structure as Pe — > cxd. As the concentration 
boundary layer wraps all the way around the disk, the downstream disturbance begins 
to look like that of a Green's function source located at the rear stagnation point. If one 
defines the "wake" as the region behind the disk enclosed by a given iso-concentration 
line, c = Co, then it is easy to see that the wake tends to a finite length as Pe ^ oo, even 
though its thickness tends to zero, like l/VPe. Physically, the (dimensionless) length 
L = 0(/oPe) = 0(1) is the distance traveled in the fiow downstream in x during the time 
for diffusion across the initial wake thickness, Iq = 0(l/-\/Pe). For example, the c = 0.5 
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contour ends roughly 2.3 disk diameters downstream from the rear stagnation point, as 
Pe — > oo. 

These results have relevance for more complicated, physical situations, such as the coat- 
ing of fibers from a gas flow, where quasi-steady advection-diffusion in a two-dimensional 
potential flow is a reasonable model for the growth dynamics (K. G. Kornev, private 
communication). It is well known that the Peclet number, based on the fiber diame- 
ter, controls the total growth rate, as well as the uniformity of the coating's thickness, 
for a single fiber. For Pe < 0.1, the flux to the disk (or fiber cross-section), cr(6',Pe), 
is nearly uniform, apart from a minor asymmetry due to the flow, which is accurately 
described by 15.24(1 . For Pe > 0.01, the high-Pe expansion (|4.16() shows how the flux 
profile becomes increasingly asymmetric and approaches the advection-dominated limit, 
(7 ^ 2-\/Pe/7r sin 0/2 as Pe — > oo. Our work provides an analytical description of how this 
crossover occurs for the flux profile and the total flux Ht).4|l . which may be useful 

in future analytical or numerical studies of the coating process. 

By quantifying the crossover in the concentration field, we also provide some insight 
into the possible effect of interactions between multiple nearby fibers during real coating 
process. In the "cloud" regime at low Pe, the concentration field approximately satisfies 
Laplace's equation, which means that other fibers in all directions can strongly influence 
the local flux, due to the long-range decay of the concentration. In the 'wake' regime 
at high Pe, the concentration remains uniform (outside a thin boundary layer) in all 
directions except directly downstream, where a thin wake forms of 0(l/-\/Pe) thickness 
and 0(1) extent (at the scale of a few particle diameters). When the mean fiber spacing 
is larger than the typical wake length, there are negligible interactions, but, whenever a 
fiber ends up in the wake of another fiber, its coating becomes much thinner in a localized 
region, which can be undesirable. Our analysis shows that this important crossover occurs 
at a critical value of Pe « 60. 

Another physical application of our results is to advection-diffusion-limited aggregation 
l(Bazant et alX 2003), the stochastic analog of the fiber coating problem, which involves 
discrete, diffusing particles collecting on a seed in a potential flow. The ADLA model may 
have relevance for the dendritic growth of mineral oxide deposits in very thin cracks in 
rocks, where Hele-Shaw flow, laden with oxygen diffusing from the surface, is believed to 
occur over geological time scales (G. R. Rossman, private communication). To describe 
ADLA growth from a finite seed in a uniform background flow, the probability measure 
for adding particles on a growing aggregate in the z plane is given by the flux to the unit 
disk for the same problem in the w plane, 

p(z, t)\dz\ = p{e, t)de cx (J^{e; Pe(i)), (7.2) 

where the renormalized Peclet number, Pe(t) = ^i(t)Peo, is given by the time-dependent 
conformal radius. Our accurate numerical solution for a^iO^Ve) has already been used 
as the growth measure in the simulations of IBaz ant et al.\ (2003). Future simulations 
or mathematical analysis could start from our nearly exact formula, cr'^'^°""^ from (|6.1|) . 
which captures the full crossover between diffusion-dominated and advection-dominated 
dynamical fixed points with increasing Pe(i). 
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Appendix A. Invertibility of the Kernel 

In this appendix we show that the symmetrized kernel, i^o(Pc|a; — a:'|), of the integral 
equation (|2.9|l is positive definite. Using the relation 

/ N 1 /"°° , cos(xt) , . , , 



the kernel is recast to the form 



cos(<;^) cos(';o^) + sin(<;^) sin(^ot) 



It follows that, for /3 > 0, 

d,of^ d, r dt cos(^ot) + sin(a) sin(^ot) ^(.^^^(.^ ^ r gjtf + h{tf ^ 

-f3 J-p J-oo \/l + J-oc \/l + 

(A3) 

by which the kernel is positive definite. The functions g(t) and h{t) are real, and given 

by 

g{t) + ih{t)^ f{u)e''^'du. (A 4) 
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